set.seed(51)
# load(file = '~/Data/Medicare/Linked AQS-Medicare/HEI PM10 Accountability/Original HEI Report Submission/PM10_withps1990.Rda')
load(file="~/PM10/pm10_withps1990_heiorig.Rda")
######## DO NOT Prune observations that are out of range ###########
# dat=subset(dat, outofrange==0)
## Redefine pscat with unpruned observations
nblock   <- 5
breaks  <- as.numeric(quantile(dat$ps, probs = seq(from=0, to=1, by = 1/nblock)))
dat$pscat  <- as.numeric(cut(dat$ps, breaks= breaks, include.lowest=T))
dat$pscat = as.factor(dat$pscat)


######## Set the Working Directory so that output is saved in the right place ######## 
wd <- "~/PM10/20150904_HEIorig_noprune_pollution/"
setwd(wd)

##############################################################
################ Specify the Array Parameters ################
##############################################################
ytrans = c("log", "change", "95") 

##############################################################
###################### Array Starts Here #####################
##############################################################
j = as.numeric(Sys.getenv('SLURM_ARRAY_TASK_ID'))
print(j)


library(splines)
library(HEIfunctions)
library(coda)
library(corpcor)

######## Choose how to transform pollution ######## 
if (ytrans[j]=="log")
  dat$y = log(dat$pm10fu)   
if (ytrans[j]=="change")
  dat$y = dat$pm10fu - dat$pm10base1990
if (ytrans[j]=="95")
  dat$y = log(dat$pm10fu95)



######## Choose how to specify the model ######## 
### Big Models ###
vars = c("Monitor", "a", "Median_income", "HS_rate", "Urban_rate", "Migration_5_year_rate", "Hispanic_cens", "White_cens", "Black_cens",
         "Current_Smoking_rate_weighted", "Female_cens", "mean_age.2001", "Female_rate.2001", "White_rate.2001", "Black_rate.2001",
         "Tot_Beneficiary.2001", "pm10base1990", "housedens", "temperature", "logpop")
## Use all the variables listed above except the census ones that are duplciated in medicare
my=glm(y~pscat + pm10base1990 + Tot_Beneficiary.2001 + Median_income + HS_rate + Urban_rate + Migration_5_year_rate + Current_Smoking_rate_weighted + mean_age.2001 + Female_rate.2001 + White_rate.2001 + Black_rate.2001 + housedens + temperature + logpop, 
       data=dat)
formula = my$formula

coords 	<- cbind(dat$Longitude, dat$Latitude)

phistart 	<- makephi(coords, 10)
philower 	<- makephi(coords, 4)
phiupper 	<- makephi(coords, 30)

tuning 	<- list(A = 0.1, psi = 0.2, theta = 1)
prior 		<- list(KIG = rep(.5,2), psi = rep(.5,2), theta1 = rep(philower, 2), theta2 = rep(phiupper,2))
starting 	<- list(B=NULL, A = NULL, psi = NULL, theta = phistart)

mod 		<- mvpsmod(formula, dat, "a", coords, nsamp=50000, thin=1,
			tuning=tuning, prior=prior, starting=starting, 
      outputbinsize = 500, outputfilename = paste("pollutionmodel_temp_", ytrans[j],".Rda", sep=""))

save(mod, file = paste("pollutionmod_", ytrans[j], ".Rda", sep=""))





